import cvxpy as cvx
import numpy as np
import matplotlib.pylab as plt
from typing import Tuple
import altair as alt
from vega_datasets import data
import pandas as pd
import plotly.graph_objects as go
# create problem data
N = 100
# create an increasing input signal
xtrue = np.zeros(N)
xtrue[0:40] = 0.1
xtrue[50] = 2
xtrue[69:80] = 0.15
xtrue[80] = 1
xtrue = np.cumsum(xtrue)
# pass the increasing input through a moving-average filter
# and add Gaussian noise
h = np.array([1, -0.85, 0.7, -0.3])
k = len(h)
yhat = np.convolve(h, xtrue)
# Data ported from matlab file
y = yhat[0:-3] + np.array(
[
-0.43,
-1.7,
0.13,
0.29,
-1.1,
1.2,
1.2,
-0.038,
0.33,
0.17,
-0.19,
0.73,
-0.59,
2.2,
-0.14,
0.11,
1.1,
0.059,
-0.096,
-0.83,
0.29,
-1.3,
0.71,
1.6,
-0.69,
0.86,
1.3,
-1.6,
-1.4,
0.57,
-0.4,
0.69,
0.82,
0.71,
1.3,
0.67,
1.2,
-1.2,
-0.02,
-0.16,
-1.6,
0.26,
-1.1,
1.4,
-0.81,
0.53,
0.22,
-0.92,
-2.2,
-0.059,
-1,
0.61,
0.51,
1.7,
0.59,
-0.64,
0.38,
-1,
-0.02,
-0.048,
4.3e-05,
-0.32,
1.1,
-1.9,
0.43,
0.9,
0.73,
0.58,
0.04,
0.68,
0.57,
-0.26,
-0.38,
-0.3,
-1.5,
-0.23,
0.12,
0.31,
1.4,
-0.35,
0.62,
0.8,
0.94,
-0.99,
0.21,
0.24,
-1,
-0.74,
1.1,
-0.13,
0.39,
0.088,
-0.64,
-0.56,
0.44,
-0.95,
0.78,
0.57,
-0.82,
-0.27,
]
)
xtrue
Formulate the problem of finding the maximum likelihood estimate of x, given y, taking into account the prior assumption that x is nonnegative and monotonically nondecreasing, as a convex optimization problem.
$ \begin{align} \mathbf{maximize} \; &\, \; p(x) = -(m/2) \cdot log(2 \pi \sigma^2) -\frac{1}{2\sigma^2} \biggl\lvert\biggl\lvert \; \sum_{\tau=1}^k {h(\tau) x(t-\tau)} -y \;\biggr\rvert\biggr\rvert_2 \newline \mathbf{subject \; to}: \, & \; x_i>=0 \newline & x_{i+1} > x_i \end{align} $
This is equivalent to:
$ \begin{align} \mathbf{minimize} \; &\, \; \hat{y} -y \newline \mathbf{subject \; to}: \, & \; x_i>=0 \newline & x_{i+1} > x_i \end{align} $
where $\hat{y} = \biggl\lvert\biggl\lvert \; \sum_{\tau=1}^k {h(\tau) x(t-\tau)}\;\biggr\rvert\biggr\rvert_2 $
# Define variables
x_constrained = cvx.Variable(N)
x_free = cvx.Variable(N)
# Define objectives
obj_constrained = cvx.Minimize(
cvx.norm(cvx.conv(h, x_constrained)[0:-3] - y[:, np.newaxis])
)
obj_free = cvx.Minimize(cvx.norm(cvx.conv(h, x_free)[0:-3] - y[:, np.newaxis]))
# Define constraints
constraints = [
x_constrained >= 0,
x_constrained[1:] >= x_constrained[:-1]
]
# Form and solve the constrained problem
prob_constrained = cvx.Problem(obj_constrained, constraints)
prob_constrained.solve() # Returns the optimal value.
# Form and solve the free problem
prob_free = cvx.Problem(obj_free)
prob_free.solve() # Returns the optimal value.
print("optimal x_constrained", x_constrained.value)
print("optimal x_free", x_free.value)
# Plot results
df_true = pd.DataFrame({"t": range(100), "x": xtrue, "type": "true"})
df_constrained = pd.DataFrame(
{"t": range(100), "x": x_constrained.value, "type": "constrained"}
)
df_free = pd.DataFrame({"t": range(100), "x": x_free.value, "type": "free"})
df = pd.concat([df_true, df_constrained, df_free])
alt.Chart(df).mark_line().encode(
x="t:T",
y="x:Q",
color=alt.Color(
"type",
scale=alt.Scale(
domain=["true", "constrained", "free"], range=["red", "green", "blue"]
),
),
).properties(width=800, height=450)
Two investments are made, with random returns $R_1$ and $R_2$. Find maximum $p_\mathrm{loss} = prob(R_1 + R_2 <= 0)$.
$R_1$ and $R_2$ have Gaussian marginal distributions with known means $\mu_1$ and $\mu_2$ and known standard deviations $\sigma_1$ and $\sigma_2$.
$R_1$ and $R_2$ are correlated with correllation coefficient $\rho$:
\begin{align} E(R_1 - \mu_1)(R_2 - \mu_2) = \rho \sigma_1 \sigma_2 \end{align}Check joint Gaussian distribution first:
mu_1 = 8
mu_2 = 20
sigma_1 = 6
sigma_2 = 17.5
rho = -0.25
from scipy.stats import norm
p_loss = norm.cdf(
0,
mu_1 + mu_2,
np.sqrt(sigma_1 * sigma_1 + sigma_2 * sigma_2 + 2 * rho * sigma_1 * sigma_2),
)
p_loss
Now: Do not assume a joint gaussian distribution. The joint distribution can be anything as long as the marginal distributions are Gaussian:
$$ \begin{align} \mathbf{maximize} \; &\, \; \sum_{R_1 + R_2 <= 0} P_{ij} \\ \mathbf{subject \; to}: \, & \; P_{ij}>=0 \;\;\;\;\;\; \mathrm{(this\ is\ a\ pdf)} \\ & \mathbf(P) * \mathbf{1} = p^{(1)} \;\;\;\;\;\; \mathrm{(marginal\ distribution\ for\ p^{(1)}\ is\ Gaussian)} \\ & \mathbf(P)^T * \mathbf{1} = p^{(2)} \;\;\;\;\;\; \mathrm{(marginal\ distribution\ for\ p^{(2)}\ is\ Gaussian)} \\ & (r - \mu_1 \mathbf{1})^T P (r - \mu_2\mathbf{1}) = \rho \sigma_1 \sigma_2 \;\;\;\;\;\; \mathrm{(correlation)} \end{align} $$Assuming a discretization of $R_1$ and $R_2$ between $r_1 = -30$ and $r_n = 70$, the marginal distributions for $p^{(1)}$ and $p^{(2)}$ are given by:
$$ \begin{align} p_i^{(1)} = \mathrm{prob}(R_1 = r_i) &=& \frac{\exp( -(r_i-\mu_1)^2 / (2\sigma_1^2))} {\sum_{j=1}^n\exp(-(r_j-\mu_1)^2/(2\sigma_1^2))}\\\,\\ p_i^{(1)} = \mathrm{prob}(R_1 = r_i) &=& \frac{\exp( -(r_i-\mu_2)^2 / (2\sigma_2^2))} {\sum_{j=1}^n\exp(-(r_j-\mu_2)^2/(2\sigma_2^2))} \end{align} $$# Compute linear grid for approximate solution
m = 100
r = np.linspace(-30, 70, m)
# Marginal probabilities:
p_1 = np.exp(-(r - mu_1) * (r - mu_1) / 2 / sigma_1 / sigma_1) / np.sum(
np.exp(-(r - mu_1) * (r - mu_1) / 2 / sigma_1 / sigma_1)
)
p_2 = np.exp(-(r - mu_2) * (r - mu_2) / 2 / sigma_2 / sigma_2) / np.sum(
np.exp(-(r - mu_2) * (r - mu_2) / 2 / sigma_2 / sigma_2)
)
# Define loss mask as indicator matrix for entries that denote loss situations (i.e. r1 + r2 <= 0)
loss_mask = np.sum(np.meshgrid(r, r), axis=0) <= 0
# Define cvx problem
P = cvx.Variable((m, m))
objective = cvx.Maximize(cvx.sum(cvx.multiply(P, loss_mask)))
constraints = [
P >= 0,
P @ np.ones(m) == p_1,
P.T @ np.ones(m) == p_2,
(r - mu_1) @ P @ (r - mu_2) == rho * sigma_1 * sigma_2,
]
problem = cvx.Problem(objective, constraints)
# Solve problem
problem.solve()
# Define points for plots
x = np.repeat(r, m)
y = np.tile(r, m)
z = P.value.flatten()
# Create 3D Mesh plot
fig = go.Figure(data=[go.Mesh3d(x=x, y=y, z=z, color="red", opacity=0.50)])
fig.show()
# Create contour plot
fig = go.Figure(data=go.Contour(x=x, y=y, z=z))
fig.update_layout(
autosize=False,
width=800,
height=800,
)
fig.show()